
# Figure 2
sims <- 10000
lm1_coefs <- mvrnorm(sims, lm1$coefficients, vcovHC(lm1, type="HC4"))
lm1_treatB <- lm1_coefs[,2]
lm1_treatC <- lm1_coefs[,3]
pe1 <- c(mean(lm1_treatB), mean(lm1_treatC))
lwr1 <- c(quantile(lm1_treatB, probs=0.025), quantile(lm1_treatC, probs=0.025))
upr1 <- c(quantile(lm1_treatB, probs=0.975), quantile(lm1_treatC, probs=0.975))

lm2_coefs <- mvrnorm(sims, lm2$coefficients, vcovHC(lm2, type="HC4"))
lm2_treatB <- lm2_coefs[,2]
lm2_treatC <- lm2_coefs[,3]
pe2 <- c(mean(lm2_treatB), mean(lm2_treatC))
lwr2 <- c(quantile(lm2_treatB, probs=0.025), quantile(lm2_treatC, probs=0.025))
upr2 <- c(quantile(lm2_treatB, probs=0.975), quantile(lm2_treatC, probs=0.975))

lm3_coefs <- mvrnorm(sims, lm3$coefficients, vcovHC(lm3, type="HC4"))
lm3_treatB <- lm3_coefs[,2]
lm3_treatC <- lm3_coefs[,3]
pe3 <- c(mean(lm3_treatB), mean(lm3_treatC))
lwr3 <- c(quantile(lm3_treatB, probs=0.025), quantile(lm3_treatC, probs=0.025))
upr3 <- c(quantile(lm3_treatB, probs=0.975), quantile(lm3_treatC, probs=0.975))

trace1 <- ropeladder(x=pe1, lower=lwr1, upper=upr1,
                     sublabels=c("Treat 1 effect", "Treat 2 effect"), 
                     sublabelsyoffset="-0.05",
                     size=0.5, lex=1.5, lineend="square",
                     plot=1)

trace2 <- ropeladder(x=pe2, lower=lwr2, upper=upr2, 
                     sublabels=c("Treat 1 effect", "Treat 2 effect"), 
                     sublabelsyoffset="-0.05",
                     size=0.5, lex=1.5, lineend="square",
                     plot=2)

trace3 <- ropeladder(x=pe3, lower=lwr3, upper=upr3, 
                     sublabels=c("Treat 1 effect", "Treat 2 effect"), 
                     sublabelsyoffset="-0.05",
                     size=0.5, lex=1.5, lineend="square",
                     plot=3)

sigMark1 <- pe1
is.na(sigMark1) <- (lwr1>0 | upr1<0)
traceSig1 <- ropeladder(x=sigMark1, col="white", group=1, plot=1)

sigMark2 <- pe2
is.na(sigMark2) <- (lwr2>0 | upr2<0)
traceSig2 <- ropeladder(x=sigMark2, col="white", group=1, plot=2)

sigMark3 <- pe3
is.na(sigMark3) <- (lwr3>0 | upr3<0)
traceSig3 <- ropeladder(x=sigMark3, col="white", group=1, plot=3)

vertmark1 <- linesTile(x=c(0,0), y=c(0,1), plot=1)
vertmark2 <- linesTile(x=c(0,0), y=c(0,1), plot=2)
vertmark3 <- linesTile(x=c(0,0), y=c(0,1), plot=3)

tile(trace1, trace2, trace3, traceSig1, traceSig2, traceSig3,
     vertmark1, vertmark2, vertmark3,
     RxC = c(3,1),
     limits=c(-33,33), gridlines=list(type="xt"),
     topaxis=list(add=T, at=c(-30,-20,-10,0,10,20,30)),
     plottitle=list(labels=c("Full Sample", "Wando County", "Yeoju City")),
     xaxistitle=list(labels="Effect size"),
     width=list(spacer=3),
     height=list(plottitle=2, xaxistitle=3, spacer=5),
     output=list(outfile="figure3", width=10))

